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ABSTRACT 

Porosity evolution of dust aggregates is crucial in understanding dust evolution in protoplanetary disks. In 
this study, we present useful tools to study the coagulation and porosity evolution of dust aggregates. First, 
we present a new numerical method for simulating dust coagulation and porosity evolution as an extension of 
the conventional Smoluchowski equation. This method follows the evolution of the mean porosity for each ag- 
gregate mass simultaneously with the evolution of the mass distribution function. This method reproduces the 
results of previous Monte Carlo simulations with much less computational expense. Second, we propose a new 
collision model for porous dust aggregates on the basis of our A^-body experiments on aggregate collisions. As 
the first step, we focus on "hit-and-stick" collisions, which involve neither compression nor fragmentation of 
aggregates. We first obtain empirical data on porosity changes between the classical limits of ballistic cluster- 
cluster and particle-cluster aggregation. Using the data, we construct a recipe for the porosity change due to 
general hit-and-stick collisions as well as formulae for the aerodynamical and collisional cross sections. Our 
collision model is thus more realistic than a previous model of Ormel et al. based on the classical aggregation 
limits only. Simple coagulation simulations using the extended Smoluchowski method show that our collision 
model explains the fractal dimensions of porous aggregates observed in a full A^-body simulation and a labo- 
ratory experiment. By contrast, similar simulations using the collision model of Ormel et al. result in much 
less porous aggregates, meaning that this model underestimates the porosity increase upon unequal-sized col- 
lisions. Besides, we discover that aggregates at the high-mass end of the distribution can have a considerably 
small aerodynamical cross section per unit mass compared with aggregates of lower masses. This occurs when 
aggregates drift under uniform acceleration (e.g., gravity) and their collision is induced by the difference in their 
terminal velocities. We point out an important implication of this discovery for dust growth in protoplanetary 
disks. 

Subject headings: dust, extinction — planetary systems: formation — planetary systems: protoplanetary disks 
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1. INTRODUCTION 

It is a widely accepted idea that the first step of planet for- 
mation in protoplanetary disks involves the collisional growth 
of submicron/micron dust grains into macroscopic aggre- 
gates. A standard scenario is that dust grains coagulate by mu- 
tual sticking, gradually settle towar d the midplane of the disk, 
and f orm a dense dust layer (e.g., IWeid enschilling & Cuzzil 
1993). It is still an open issue whether subsequent plan- 
etesimal formatio n is achieved by the gravitational instabil- 
ity of the layer dSafronovl 11969): iGoldreich & Ward! 119731: 
ISekivalfT99 8: Youdin & Shull2002l) or by the direct coll i sional 
growth of the aggreg ates (WeidenschiHin g~& Cuzzil 119931: 
Stepin ski & Valageaslll997t [Brauer et al.l 120081) . To address 
this issue, further understanding on earlier evolutionary stages 
is needed. 

Most of the previous studies have simplified dust aggregates 
as compact, nonporous spheres. However, both numerical and 
laboratory experiments have revealed that aggregates are not 
at all compact, but have an open, flu ffy structure (for a review, 
see Blum 2004; Dominik et al. 2007). This is particularly true 
for aggregates formed at an early growth stage where the col- 
lisional velocity is so low that collision al compression is n eg- 
ligible. It has been ob served in jV-bodv (iKempf et al.ll9 99) as 
well a s experimental dWurm & Blumll l998l:lBlum etal.ll 19981 
2000) studies that the outcome is an ensemble of fractal aggre- 
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gates with the fractal dimension of D < 2. This fractal growth 
lasts until the impact en ergy becomes high e nough to cause 
collisional compaction (Suva maet al.| [2008). The porosity 
change due to compressive as well as destructive collisions 
are al s o studied by so me aut hors (e.g.JWada et all2007l l2008l 
120091: ISuyamaet ail l2008t iPaszun & Dominikl 120091) using 
jV-body simulations inclu ding monomer surface interactions 
dDominik & Tielensll 1997b . 

There are several strong reasons why the porosity evolution 
of aggregates is important in studying dust evolution in proto- 
planetary disks. These include the following. 

1 . The porosity affects the aerodynamical property of ag- 
gregates. In a gas disk, the motion of a small aggre- 
gate is controlled by the friction force from the ambient 
gas. The friction coefficient can vary by many orders of 
magnitude depending on whether the aggregate is com- 
pact or fluffy. For example, a fractal (D < 2) aggregate 
made of one thousand grains receives a friction force an 
order of magnitude stronger than a compact aggregate 
of the same mass. This difference is significant when 
one considers the formation of sedimentary dust layer 
on the disk midplane. 

2. The porosity even affects the turbulence of proto- 
planetary disks. Magne torotational instability (MRI; 
Balbus & Hawleyl 119911) . which is recognized as the 
most likely mechanism driving disk turbulence, oper- 
ates only in a re gion where the ionization degree is 
sufficiently high (Gammie 1996). The ionization de- 
gree strongly depends on the surface area of dust ag- 
gregates, since the recombination of ionized gases ef- 
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ficiently takes place on dust surfaces. Again, a frac- 
tal (D < 2) aggregate made of one thousand grains 
absorbs ionized gases an order magnitude more effi- 
ciently than a compact grain of the same mass. For 
this reason, compact dust growth tends t o make the gas 
environment turbulent (Sano et al.l 120001) . while frac- 
tal dust growth tend s to keep the environment lami- 
nar (Okuzumi 2009). This difference is critical to the 
evolution of dust itself, since the turbulence causes the 
diffusion and collisional fragmentation of aggregates 
dBrauer et al.ll2008tlJohansen et al.ll2008l) . 

3. The fractal nature of aggregates matters when we try to 
interpret observed emission spectra of protoplanetary 
disks. For example, the 10/im spectral feature charac- 
teristic of submicron silicate grains quickly vanishes in 
compact aggre gation, but onl y slowly in fractal (D ~ 2) 
aggregation (Min et al. 2006). This means that the con- 
ventional compact dust model generally underestimates 
the degree of dust growth in the observed objects. 

Thus, the consistent modeling of the growth and porosity evo- 
lution of dust aggregates is needed to better understand and 
predict dust evolution in protoplanetary disks. 

Theoretically, however, the consistent modeling is not an 
easy task. First of all, we must consider the evolution of ag- 
gregates in the two-dimensional (2D) parameter space of mass 
and porosity. However, such 2D calculation is generally much 
more elaborate than one-dimensional (ID; i.e., mass only) 
one. Furthermore, such calculation requires a reliable colli- 
sion model for porous aggregates. Here, a "collision model" 
refers to a set of (1) the definition of the "porosity," or "vol- 
ume," of an aggregate, (2) a recipe that determines the poros- 
ity change due to collisions, and (3) formulae for collisional 
and aerodynamical cross sections that determine the collision 
rate of porous aggregate pairs. To obtain a reliable collision 
model, one needs a sufficient amount of empirical data on ag- 
gregate collisions. 

Because of the above theoretical difficulty, there are only 
a few studies addressing the coagulation of porous dust ag- 
gregates in the literature. A 2D simulation on porous dust 
growth was first done by Ossenkopf (1993) in the context 
of dust extinction in molecular clouds. He used a colli- 
sion model based on the A^-body simulations of the ballis- 
tic cluster- cluster and pa rticle-cluster aggregation (BCCA 
and BP CA; lMeakinlfT99l . Ormel et al. (2007; hereafter, 
IOST07I) have presented a Monte Carlo method for study- 
ing porosity evolution of dust aggregates in protoplanetary 
disks. In this method, porous aggregates are represented 
by the same number of computational particles with mass 
and porosity, and the collisions among the particles are suc- 
cessively calculated using random numbers. The change in 
porosity on each individual collision is determined from a 
recipe base d on previous jV-body experiments on agg regate 
collisions (Ossenkopfl ll993t iDominik & Tielens 1997). The 
Monte Car lo approach has been fu rther developed indepen- 
dently by lOrmel & Spaansl (12008b and IZsom & Dullemondl 
(2008) to achieve a high dynamical range in the parame- 
ter space. The compression/f ragmentation mod el of IOST07I 
has been recently revised by Orm el et al.l (120091) on the basis 
of numerical experiments on low-mass aggregate collisions 
dPaszun & Dominikll2009l) . 

In this study, we provide new theoretical tools useful 
for studying the collisional evolution of porous dust aggre- 
gates. First, we present a new numerical method to solve the 



growth and porosity evolution of aggregates. This method 
is an extension of the conventional Smoluchowski method, 
a method c ommonly used for ID simulations of dust aggre- 
gates (e.g-lTanaka et alJl2005b iDullemond & Dominikll2005l : 
Brau eretalJl2008h . The core of our new method is the ap- 
proximation that the porosity distribution is narrow for each 
aggregate mass. With this approximation, we derive a closed 
set of moment equations of the 2D Smoluchowski equation 
for the mass distribution function and the averaged mass- 
porosity relation. These equations are formally similar to 
the conventional ID Smoluchowski equation, and makes it 
possible to adopt well-established algorithms for ID prob- 
lems. We confirm that this approach ind eed repr oduces the 
results of the Monte Carlo simulations bv lQST07l with much 
less computational effort. Thus, our method will be particu- 
larly useful for adding the porosity evolution to glob al simu- 
lations including radial drift (e.g., Brau er et al.ll2008l) or cou- 
pled simulations inv olving hydrodynamic calculation (e.g., 
iJohansen et al J 120081) . Second, we provide a new collision 
model based on our A^-body experiments for various types 
of aggregate collisions. As the first step of the modeling, 
we focus on "hit-and-stick" collision, i.e., low-velocity col- 
lision involving neither compression nor fragmentation. In a 
typical protoplanetary disk, the hit-and-stick picture well rep- 
rese nts the growth of du st aggregates up to several centime- 
ters dSuvama et al. 2008). In contrast to previous studies, we 
first use empirical data on the porosity evolution between the 
BCCA and BPCA limits. We will see that the numerical sim- 
ulations of porosity evolution using our collision model well 
agree with the results of previous full A^-body and laboratory 
experiments. Extension of our model to compressive and de- 
structive collisions will be made in the future work. 

This paper is organized as follows. In Section 2, we de- 
scribe our extended Smoluchowski method and its numerical 
implementation. In Section 3, we use this method to tr y to re- 
produce the results of the Monte Carlo simulations by OST07. 
In Section 4, we present a new porosity model as well as the 
results of our A^-body simulations from which we construct 
the porosity increase recipe. In Se ction 5, we compare this 
collision model with that of OST07. Section 6 is devoted to 
the summary. 

2. THE EXTENDED SMOLUCHOWSKI METHOD FOR 
POROUS DUST GROWTH 

In this section, we describe a new method to solve the co- 
agulation and porosity evolution of dust aggregates. A com- 
parison between our method and the Monte Carlo method is 
made in Section 3. 

2. 1 . The Multi-dimensional Smoluchowski Equation 

We denote the internal structure of each aggregate by a set 
of parameters I = {I m ,I <2 \. . .}. Conventional coagulation 
models adopted the mass M of an aggregate as a unique pa- 
rameter, i.e., I = {M}. This study treats the volume V as an 
additional parameter for porous aggregates. We assume that 
the internal structure of an aggregate created by a collision 
is uniquely determined by those of the collided ones. This 
means that the structure parameter of the new aggregate, Ij+2 
is a function of those of the old ones, Ii and I2. 

We denote the distribution function of aggregates by /(I). 
The evolution of /(I) is determined by the multidimensional 
Smoluchowski equation 

^M= l -Jd\'d\"K{t;\")f{V)f{\") 
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x5[I 1+2 (I';I")-I] 
-/(I) j dl' K&IW), 
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where K(Ji',lj) = K(If,Ii) is the collision kernel defined as the 
product of the collisional cross section cr co u and the relative 
speed Am of colliding pairs. The first and second terms in 
Equation ([1) represent the "gain" and "loss" of aggregates 
with I due to collisions, respectively. 

If I only contains M, Equation {T|i reduces to the conven- 
tional, ID Smoluchowski equation 



dfm = i 

dt 2 



/ dM' K(M';M-M')f(M')f(M-M') 
Jo 

POO 

-f(M) / dM' K(M;M')f(M'), (2) 
Jo 

where we have used the mass conservation M\ +2 (M' ,M") = 
M' + M". It is easy to check that Equation (O ensures the 
conservation of the total mass density p c i = J Mf(M)dM. 

Now we consider the 2D case, I = (M, V). In this case, equa- 
tion (]]) is written as 



9/m,v 
dt 



dM' 



dV'dV" K^mi yn /m' .v'/m-m 1 y" 



x5[(v 1+2 )^: v „ 



-v\ 



-fu 



■ J dM'dV' K™; v v ,f M , r 



(3) 



where f u ,v = f(M,V), K%^„ = K(M> ,V';M" ,V"), and 

(Vi+2)m»P = V 1+ 2(M',V';M",V"). Again, one can easily 
check that Equation © ensures the conservation of the total 
mass density p d = j Mf(M,V)dMdV for arbitrary V\ +2 . 

2.2. The Volume-averaging Approximation 

In principle, the 2D Smoluchowski equation ([3]l can be 
solved once K and V\+ 2 are given. In practice, however, direct 
numerical integration of Equation (01 requires a huge com- 
putational expense. To see this, let us consider that we try 
to solve Equation (f3]i by dividing each parameter dimension 
with N fixed bins. At each step, we have to calculate colli- 
sion events over all pairs of the bins. For a P-dimensional 
problems, the number of the fixed bins is M v , so the total 
number of the collision pairs is approximately Af 2V . Thus, 
we see that a 2D calculation of the Smoluchowski equation is 
approximately N 2 times heavier than an ID one. Of course, 
one is free to take very small M, but then it would result 
in a very poor resolution of f(M,V). For a practical use, 
some prescription is clearly needed to manage both a rea- 
sonably small expense and reasonably good accuracy. One 
option may be to use fewer bins and instead continuously 
adjust them so that they span a parameter region with large 
f{M,V). This approach appears to be efficient only if the ad- 
justment requires only a few calculations. Anoth er optio n is 
to adopt a direct Monte Ca rlo method, as done by OST071and 
IZsom & Dullemondl (120081) . Monte Carlo methods are effec- 
tively similar to the moving bin method since the Monte Carlo 
particles are continuously redistributed in the 2D space. 

Now we propose a new method. A basic strategy of this 
method is to give up following the details of the volume dis- 
tribution and to only follow its average for each mass. This 



enables us to calculate the porosity evolution with only 0(Af 2 ) 
calculations at each time step as is for ID problems. 

For a rigorous formulation, we introduce the moments of 
the 2D distribution function 



n(M)= / f(M,V)dV, 



V(M) = 



n(M) 



Vf{M,V)dV, 



(4) 



(5) 



where n(M) and V(M) denote the number density and mean 
volume of aggregates with mass M, respectively. The evolu- 
tion of these quantities is described by the moment equations 
of the 2D Smoluchowski equation, Equation Taking the 
zeroth- and first-order moments of equation ([3]) with respect 
to V, we obtain 



dn(M) _ 1 
dt ~2 



dM' K(M';M-M')n(M')n(M-M') 



/"OC 

-n(M) / dM'K(M;M')n(M'), (6) 
Jo 



and 



d[V(M)n(M)] 1 



<:)t 



- / dM'Vi +2 K(M';M-M') 
2 Jo 

xn(M')n(M-M') 

-n(M) / dM' VK(M;M')n(M'), (7) 
Jo 



where K, Vi +2 K, and VK are the volume averages of K, Vi +2 K , 
and VK defined by 

K(M';M") f dV'dV" <-;: / ^' / ^" (8) 
J m > v n(M') n(M") 

V^K(M';M")^ J dV'dV" (V^^K^,, 



, fwy fM",v" 
' n(M') n(M") ' 



(9) 



VK(M;M') = I dVdV VE£\, (10) 
J n(M) n(M') 

Note that the above volume-averaged quantities generally 
depend on higher-order moments of V. To solve Equations 
(O and ^} in a closed way, some additional assumption about 
the higher-order moments is needed. 

Now we try to close Equations (|6]l and Q in terms of n(M) 
and V(M) by making a simple assumption. The simplest one 
is that the volume distribution is so narrow that the volume 
moments higher than the first order can be ignored. This as- 
sumption is equivalent to supposing that the distribution func- 
tion f(M,V) is well approximated in the form 



f(M,V) = n(M)5[V-V(M)l 



(11) 



Using this assumption, the volume integral in equations ([8])- 
( fTOb can be performed for arbitrary K to give 

K(M';M") = K(M' ,V(M');M" ,V(M")) , (12) 



Vi +2 K(M';M") = Vi +2 (M' ,V(M');M" ,V(M"))K(M';M"), 

(13) 
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and 



VK(M;M') = V(M)K(M';M"). 



(14) 



Substituting these into equations (JSJ and ((7J, we obtain the 
closed forms of the moment equations 



dn(M) _ 1 
dt ~2 



M 



dM K 



I „M',V(M') 



M-M' .y {M-M' ) 



,n(M')n(M-M') 



(15) 



and 



9lV(M)n(M)] = l_ ^ M , V(M/) 

g t 2 J [+z 'm-m'.v(m-m') 

* K ul%L M MM')n{M-M') 

M-M' ,V (M-M' ) ' v J 

poo _ 

-V(M)n(M)J dM' K^ f) n(M'). 

(16) 

In this way, we have reduced the 2D Smoluchowski equation, 
Equation ||3}, to a closed set of two ID equations (T5[ and ( fTol l 
just by imposing a simple approximation, Equation ( fTTb . We 
refer to this assumption as the volume-averaging approxima- 
tion. 

The set of Equations (TT3T > and (fT6b has several advantages 
over the original full 2D equation (0. First, since all the vol- 
ume integrals has been already performed, Equations ( fT3T > and 
( fT6l ) only requires (D(Af 2 ) calculations at each time step. Sec- 
ond, as explained in the following subsection, one can solve 
these equations in just the same way as one solves the conven- 
tional ID equation The drawback is that we do not know 
in advance whether the approximation made by Equation ( fTTT i 
is reasonable or not for a coagulation problem to be solved. 
Nevertheless, we will see in Section 3 that this approximation 
is surprisingly successful in reproducing the results of previ- 
ous full 2D calculations. 

It is worth mentioning here that the above formulation can 
be extended to other aggregate parameters. For example, 
it is straightforward to derive equations like Equations $15[ 
and (O with the electric charge Q as an additional param- 
eter. This charge-averaging approximation will be particu- 
larly useful in a situation where aggregate collisions lead to 
charge separation betwee n small and large aggregates. Re- 
cently, iMuranushil d2009l) has shown that, with a sufficiently 
high dust-to-gas ratio, the charge separation could occur for 
icy aggregates in protoplanetary disks. The charge separation 
is potentially important since it can lead to the lightning in the 
disks, which is a poss ible mechanism for chondrule formation 
dDesch & Cuzzill2~000l) . The charge-averaging approximation 
will offer a powerful tool for further investigation of this issue. 

2.3. Numerical Implementation 

Before proceeding to the numerical implementation of 
Equations ( TT3T > and fl6l . we briefly review how the ID Smolu- 
chowski equation, Equation (f2j), has been solved numerically. 
Let us rewrite Equation (f2| as 



d{Mf M ) 
dt 



M/2 



dM' M l+ iK$, fa.fi 



dM'MRMfafa,, 



(17) 



where f M = f(M),K%, =K(M';M"),M" =M-M', andM 1+2 = 
M'+M" =M. This equation means that the mass density 
Mfa increases by the collision between aggregates M' and 
M"(> M') at a rate M 1+2 K^',fa>fa", and decreases by the 
collision between M and M 1 at a rate MK^fafa' ■ A numer- 
ical calculation of a ID coagulationproblem often employs 
equation (TTTb instead of equation (Q, and regards Mfa as 
a fundamental quantity t o be evolved (e.g., Nakaga wa et alj 
ll98UlTanaka et afl l2005). This formulation ensures the con- 
servation of the total mass density p,/ = J Q Mn(M)dM in an 
exact way. 

Now we go back to the 2D case. Let us rewrite Equations 
( TBT l and (Tm as 



d(Mn M ) 
dt 



M/2 



dM' M 1+2 K^, fMl n Mtl 



(18) 



d(Vn M ) 
dt 



M/2 



dM (V U 2) M ;, ^J, , V , F y "M« n M » 
dM'VK^ Mt) n M n M ,, (19) 



respectively, where n M = n(M). Note that these equations are 
very similar to the ID equation ( TP7I ). In fact, Equation ( fT8l re- 



duces to Equation (fTTIi just by rewriting K M ,/^ M , r> as M" . 
Furthermore, Equations ( fT8l and ( fT9l ) are formally identical 
except that M and Mi+2 in the former equation are replaced 
by V and V\+ 2 in the latter. Thus, we need not to prepare any 
special numerical scheme to solve Equations ( fT8l and ( fT9l ); 
we only need to adopt a scheme that has been developed to 
solve the conventional equation (Tit . 

In this study, we solve equatio ns ( fTSI ) and (T% u sing a 
fixed-bin scheme dNakagawa et al Jl 1 98 lb lTanakaetalj r20Q5: 
iBrauer et alj|2008f) . In this scheme, the mass space is divided 
into discrete bins and collision events are calculated among 
the bins. For monodisperse monomers with mass mo, the mass 
regions mo<M < Mbd^o are divided into linearly spaced bins 
with representative mass Mk = kmo (lc= 1,2,... ,J\fbd), and the 
region M > A/i^mo are divided into logarithmically-spaced 
bins with M k = \Q k l Ml "'M k -i{k > JV bd ). For polydisperse 
monomers, the mass space is just divided logarithmically. In 
this study, we do not consider poly disperse monomers. As 
pointed out by lOhtsuki et alJ (1990), the choice of Mbd must 
be carefully done because a small value of Mbd can cause arti- 
ficial acceleration in coagulation at the high-mass distribution 
tail. iLed (|2000) performed a couple of simulations with dif- 
ferent Mbd and found that the numerical solutions converge 
for Mbd > 40. In this study, we take Mbd = 40 or 80 (or equiv- 
alently, Mk+x/Mk = 1 .06 or 1 .03) depending on problems con- 
sidered. 

Transfer of mass and volume in the mass space are cal- 
culated in the following way. The number Q,j of collision 
events per unit time between two bins i and j{> i) is writ- 
ten as Qij = Kjjn/nj for j > i and = (l/2)Kunf for j = i. 
The total mass densities transferred from bins i and j in the 
events are M/Q/j and MjQjj, and the total volume densities 
transferred are V,2o and V jQij, respectively. As a natural 
consequence of the logarithmic binning, the mass of the re- 
sulting aggregates, Mn = Mj + Mj, does not necessarily coin- 
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cide with any of the bin masses, and generally goes between 
Me and Mm with some £(> j). Therefore, we have to deter- 
mine the mass densities transferred to bins I and £+ 1 in some 
way. I n this study, we adopt a prescription o flKovetz & Olundl 
( 1969), in which the mass densities transferred to bins I and 
£+1 are determined as eM/jQ/j and (l-s)MijQij respectively, 
where e = {Mi+\ - Mij)/(Mi+\ - Mi). This prescription en- 
sures the conservation of the total number and t otal mass of 
the re sulting aggregates (see Appendix A. 1 of Brau er etail 
2008). Similarly, the volume densities transferred to bins £ 
and £+1 are determined as eVyQy and (1 — £)Vf+/<2y, where 
Vij = Vna(Mi,V(Miy,Mj,V(Mj)). For £ = j, we compute the 
gain and loss of the mass and volume densities in bin j simul- 
taneously, i.e., only compute the net gains eMjQjj and (Vy — 
V(Mj))Qjj. This simultaneous computation avoids systemati- 
cally increasing errors in the densities due to the near cancel- 
lation of the gain an d loss terms beyond the doubl e precision 
(see Appendix B of Dull emond & D ominik 2005 ). Also, we 
adopt the "active" bin method dLedl2000t iTanaka et al.ll2005l) 
in which collisions are not allowed for mass bins with negli- 
gibly small mass/number densities. This method greatly re- 
duces the computational expense because high-mass and low- 
mass bins are mostly empty i n early and la t er ev olutionary 
stages, respectively. Following ITanaka et al.l (120051) . we take 
the critical minimum mass density for each mass bin to be 
10~ 25 times t he total mass densit y, pd, of the system. A test 
simulation by Tanak a et al.l (12005b shows that this prescription 
hardly affects the evolution of the mass distribution function. 

For time integration, we adopt the explicit, fourth-order 
Runge-Kutta method. The time step At is continuously ad- 
justed so that the fractional decrease in mass density at any 
bin does not exceed a constant parameter 6 (typically be- 
tween 0.01 and 0.1). Although we use the explicit integra- 
tion scheme for s implicity, the use of an implicit scheme (e.g., 
iBrauer et alj|2008l) would further accelerate the calculation. 

3. COMPARISON BETWEEN THE EXTENDED 
SMOLUCHOWSKI METHOD AND FULL 2D MONTE 
CARLO METHODS 

To demonstrate the accuracy of our extended Smolu- 
chowski method, we apply this method to try to reprodu ce the 
results of full 2D Monte Carlo simulations by OST07. This 
problem has been also solved by iZsom & Dullemond! (120081) 
to compare their Monte Carlo method with that of lOST07t 

In Sections 3.1 and 3.2, we briefly summarize the colli- 
sio n model and the protoplanetary disk model assumed in 
the IOST07 problem. In Section 3.3, we show t he result s of 
our calculation and comp are them with those of OSTQ3 and 
IZsom & Dullemond! d2008l) . 

3.1. The Collision Model o AOSWA 

Here we summarize the porous aggregate model adopted 
by IOST07L An aggregate is modeled as a cluster consisting 
of N monodisperse dust grains (monomers) of mass mo and 
radius qq. The porosity of an aggregates is represented by the 
"enlargement factor" 

V A 

^P = ^, (20) 
where V* = NVo is the compact volume, and 

/A\ 3/2 

Va=V [ — ) (21) 



is the "volume" defined by the projected area A of the ag- 
gregate. Vo and Aq are the volume and the geometric cross 
section of a monomer given by Vo = (Air/tya^ and Ao = ira^, 
respectively. Below, we refer to Va as the area-equivalent vol- 
ume. A compact aggregate has if> = 1, while a porous aggre- 
gate has i/j > 1. Each aggregate is thus characterized by its 
mass M = Nmo and enlargement factor ip. 

The enlargement factor after a collision, Vi+2 = 
ipi +2 (Mi,ip\;M 2 ,ip 2 ), is determined in two different ways de- 
pending on the value of the impact energy. The impact energy 
is defined by E = (1 /2)M(Ak) 2 , where M = M\M 2 /{M X +M 2 ) 
is the reduced mass and Am is the relative velocity. If 
E is smaller than the critical restructuring energy E lestI , 
the collided aggregates just stick to each other without 
compression ("hit-and-stick" regime). If E > E Ksa , the 
collision involves the compaction of the aggregates and 
another formula is applied for ip\ +2 (compaction regime). The 
critical restructuring energy is set to the rolling energy E m u 
given by = (irao/2)F IO \i, where F m w is the rolling-friction 
force. IOST07I adopted F mU = 8.5 x 10" ^(7/14ergcm" 2 )d yn 
as suggested by a laboratory experiment dHeim et alj [1999). 

In the hit-and-s tick regi me, the formula for 1/11+2 is given by 
(Equation (15) of|OST07j) 

. M 1 ^ 1 +M 2 tp2 {M 1 tP 1 +M 2 i/) 2 \ 3Scca/2 - 1 

Vl+2 = 7Z 7Z T7-, + %dd, (22) 

Mi+M 2 \ MiV>i / 

where <5cca = 0.95 and ip^ is an additional factor only rele - 
vant to BPCA-like collisions (see Equation (19) of lOST07l) . 
Rewriting equation ( f22b in terms of Va instead of ip, we obtain 

V AM 2 = V A . 1 ( VA y*^' 2 ) ' +^ddV*,i+2, (23) 

where V*,i+2 = V*,i + V*,2 = Vo(N\ +N 2 ) is the "compact vol- 
ume" (the volume occupied by the constituent monomers) of 
the newly formed aggregate. OST07 constructe d the above 
formu la so that it agrees with the iV-body results of Oss enkopil 
d 19931) in the BCCA and BPCA limits. However, as we show 
in Section 4, the above formula is invalid between the two lim- 
its. For +N 2 < 40, the enlargem ent factor t[>\+2 is d irectly 
determined by the fitting formula of Ossenkopf ( 1993). 
In the compaction regime, ip\ +2 is given by 

K.i +2 (V'i + 2-l) = (l-/c)[V,.i(V'i-l) + V, ,2(^2-1)], (24) 

where f c =E/[(N\ +N 2 )E T0 \{\. Using the relation V A = V*tp, 
we obtain the corresponding formula for Vaa+ 2 , 

VA M2 = V A>1 +V A , 2 -f c (V AA + V A , 2 -V* A+2 ), (25) 

IOST07I constructed this formula by assuming that the rela- 
tive decrease in the porous volume Va — V* = V *(il>— 1) scales 
linearly with E. However, IWada et ail (12008b have recently 
shown that this formula does not reproduce the results of A^- 
body simulations. 

The collisional cross section cr co n is given by a co \\ = n(aA,i + 
0-a,2) 2 , where 




is the area-equivalent radius, i.e., the radius defined by the 
projected area A. 
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10" 3 10~ 2 10" 1 
compact size a* [cm] 

Fig. 1. — Evolution of the mass distribution for a T = 10 -4 porous 
growth model calculated using the extended Smoluchowski method 
described in Section 2. The numbers 1,...,7 represent t = 
10'yr,. . . , 10 7 yr. This figure corresponds to Figure 10C in OST07. 

3.2. Protoplanetary Disk Model 

|QST07| adopted the minimum-mass solar nebula (MMSN) 
model lHavashilll98l in which the gas surface density T, g 
and the temperature T are given by power laws T, g oc R~ 3 ^ 2 
and T oc IT 1 ' 2 , where R is the distance from the central star. 
They considered two disk positions R = 1AU and 5AU, and 
assumed rocky dust for R = 1AU and icy dust for R = 5AU. 
In this study, we consider the case of R = 1AU, and take the 
dust parameters as p g J pd = 240, arj = 0.1 /mi, po = 3.0 g cm -3 , 
and 7 = 25 erg cm" 2 . The strength of disk turbu lence i s pa- 
rameterized by the so-called a-parameter a r . IOST07I con- 
sidered three cases: a r = 10" 6 , 10" 4 , and 10" 2 . The vertical 
position is taken as one scale height above the midplane. The 
collision velocity is induced by Brownian motion and turbu- 
lence. An aggregate is removed from the population (or "rains 
out" to the midplane) if its stopping time ry exceeds the crit- 
ical value r la i n = &t/Q where ft is the Kepler rotational fre- 
quency. The mass of the central star is taken as 1M Q , so that ft 
is given by Q = 2.0 x 10" 7 rad s" 1 at R = 1 AU. IOST07I adopted 
p g = 8.5 x 10" 10 g cm" 3 , c s = 1.6 x 10 5 cm s" 1 , 3 p g /p d = 240, 
the dust material density po = 3.0 gem" 3 , and the critical 
rolling energy £,011 = 2.4 x 10~ 9 erg. 

3.3. Results and Comparison 

Following [QST07L we have performed simulations for four 
growth models: three porous (a T = 10" 6 , 10" 4 , 10" 2 ) and one 
compact (aj = 10~ 4 ). We computed the evolution of n(M) 
and Va(M) using Equations ( TT8l and ( fT9b with the numerical 
scheme described in Section 2.3. The control parameters of 
the scheme are set to Mm = 40 and 5 = 0.1. We have confirmed 
the numerical result is insensitive to the values of Mm and 

3 Note that the va lues of p s and c, adopted bv lOST07l differ from the orig- 
inal MMSN values. OST07 made some order of unity errors in calculating 
these values from £» and T. 
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a T =1QT(C) 
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10~ 5 10~ 4 10~ 3 10~ 2 10~ 1 10° 10 1 
(a,) m [cm] 

Fig. 2. — Evolution of the mean aggregate sizes for different mod- 
els, calculated using the extended Smoluchowski method. {a A ) m and 
(a*) m are the mass-weighted averages of the porous and compact 
sizes, respectively. a T denotes the turbulence parameter, and the 
letters 'C and 'P' indicate the compact and porous models, respec- 
tively. The vertical axis {a A ) m /{a*) m is essentially equal to V l/3 - The 
cross symbols indicate the average porous size a A (horizontal) and 
a A /a* (vertical) of rain-out aggregates. The numbers 1 ,. . .,7 represent 
t = 10'yr,. . . , 10 7 yr. This figure corresponds to Figure 12 in OST07. 

6 as long as Mbd > 40 and 5 < 0.1. The volume Va,\+i of 
newly formed aggregates is calculated in accordance with the 
collision model of OSTQ3 

Figure Q] shows the evolution of the mass distribution func- 
tion for the porou s, aj = 1 0" 4 model. This figure corresponds 
to Figure 10C in OST07. One can readily see that the evo- 
lution of our dis tribution function is qualitatively very simi- 
lar to that of the lOST07t nearly monodisperse growth due to 
Brownian motion for t < 10 3 yr, exponential growth driven 
by turbulence for 10 3 yr < t < 10 4 yr, and "rain-out" with 
porous size a a ~ 2 cm for t > 10 4 yr. More detailed com- 
parison shows that our distribution curves fall well inside the 
error bounds of their Mon te Carlo results. This means that our 
result agrees with that of OST07 even quantitatively. 

Figure [2] shows the evolution of the mean compact 
size (a*) and the mean porous size (a.A)m f° r differ- 
ent growth models. Here, a* = (3V*/47r) 1//3 = N^ 3 a 
is the compact (mass-equivalent) volume, and (F) m = 
J F{M)Mn{M)dM / J Mn(M)dM is the mass-weighted aver- 
age of a function F(M). The mass-weighted average approxi- 
mately corresponds to the value at the peak of the mass spec- 
trum Ma(a, t )n(fl*). The vertical axis (flA)m/ (a*) m is roughly 
equal to the third root of the mean enlarg ement fa ctor, (?/W 3 ). 
Comparing this figure with Figure 12 in lOST07l we find ex- 
ce llent agr eement between our evolution curves and those 
of OST07 except that a slight difference is seen in the late- 
time behavior for the porous, aj = 10" 2 model. The cross 
symbols in Figure [2] indicate the mean values of «a (hori- 
zontal axis) and 0^/0* (vertical axis) of the rain-out aggre- 
gates for the four gro wth models. We compare our values 
with those of OST07 in Table [TJ Here "Smol+VA" refers to 
the results of our extended Smoluchowski method with the 
volume-averaging approximation, and "MC" refers to the re- 
sults of the full 2D Monte Carlo method taken from Table 2 of 
OST07. We find that these values agree within errors of 10%. 

From the above comparison, we conclude that our ex- 
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TABLE 1 

Properties of Rain-out Particles ( 1 AU models ) 







a A (cm) 




a A la* 




Model 


Smol+VA MC 


Smol+VA 


MC 


OLT = 


1(T 4 , porous 


1.9 2.1 ±0.1 


4.4 




4.4 


OLT — 


ICT 4 , compact 


0.022 0.024 


1 




1 


OLJ = 


10~ 2 , porous 


7.6 7.8±0.7 


1.5 




1.6±0.1 


Ctj = 


10~ 6 , porous 


0.019 0.021 


4.4 




4.4 



tended Smoluchowski method well repr oduces the results of 
the Monte Carlo simulations by OST07. 



3.4. Merits and Drawbacks of the Extended Smoluchowski 

Method 

In closing this section, we point out some advantages and 
disadvantages of our extended Smoluchowski method over 
the previous Monte Carlo methods. 

The most remarkable advantage of the extended Smolu- 
chowski method is the efficiency in numerical calculations. 
The CPU time required for our method to perform each of the 
above simulations is approximately on e minu te on a 3GHz 
CPU. In contrast, IZsom & Dullemondl (12008b reported that 
their Monte Carlo method req uired 10 minutes for the same 
simulation (the calculation of IOST07I is less efficient than 
thi s since they did not use a gro uping algorithm as don e 
by IZsom & Dullemondl J2008) and lOrmel & Spaansl (T2008)). 
This means that our method is at least an order of magnitude 
more efficient than the Monte Carlo methods. The high ef- 
ficiency of our method is attributed to the volume-averaging 
approximation we introduced in Section 2.2. 

Another advantage is that we can use any numerical scheme 
having developed for the conventional Smoluchowski equa- 
tion. For example, one can further accelerate ou r method just 
using an implicit time integration developed by Brau er etall 
(2008). The implicit integration is particularly useful when 
one tries to include the fragmentation of aggregates that can 
make the problem ex tremely stiff in time dBrauer et al.ll2008l 
iBirnstiel et alJl2009h . By contrast, the implicit integration is 
inapplicable to Monte Carlo methods. 

There are two main drawbacks in our method. First, our 
method cannot be used to a problem in which the volume- 
averaging approximation is invalid. For example, the volume- 
averaging approximation will break down if a coagulation 
problem to be solved involves the fragmentation of aggregates 
and if the fragments obey broad and flat porosity distribu- 
tion for each fragment mass. For t his case, we recommend 
imp roved Monte Carlo m ethods by lOrmel & Spaansl (2008) 
and lZsom & Dullemondl (120081) . Second, our Smoluchowski 
method is inapplicable to studying "runaway" growth in 
whic h the discretenes s of the number density becomes impor- 
tant dWetherii]||1990i) . Runa way growth will be well studied 
by the Monte Carlo code of lOrmel & Spaansl (120081) . which 
can, in principle, deal with a single aggregate. 

4. A NEW HIT-AND-STICK COLLISION MODEL 

The collision model of OST07 is only based on the knowl- 
edge of the porosity evolution in the BCCA and BPCA limits. 
This means that there are no empirical supports that validate 
their model for general types of aggregation. In this section, 
we present a new collision model based on A^-body simula- 
tions for more general types of aggregation. 

As a first step toward a comprehensive model, the present 



study focuses on "hit-and-stick" collisions, i.e., collisions 
involving neither compression nor fragmentation. The hit- 
and-stick approximation is valid when the collision en- 
ergy E is smaller than the rolling-friction energy E I0 \\ 
(Domi nik & Tielensl Il997t see also Section 3.1). Assum- 
ing head-on collisions between equal-sized icy aggregates 
(made of O.l/xm monomers) at 5AU in a laminar MMSN, 
for instance, the compression and fragmentation is safely ne- 
glected until the monomer number of the aggregates becomes 
~ 10 11 , or the porous size becomes ~ cm (Su vama et alj 
2008). Modeling of porosity change due to the compression 
and fragmentation requires A^-body simulations that take into 
account the surface interaction betw een constituent particles, 
as recently done by several authors dWada et al]|2007l 120081 
I2009t ISuvama et al]|2008l: iPaszun & Dominikll2009l) At the 
present, this kind of simulations have been only performed 
for limited types of collisions (e.g., equal-sized, head-on, or 
low-mass collisions), and we thus have no sufficient empir- 
ical data on how the compression and fragmentation affects 
the porosity of resultant aggregates for more general cases. 
For this reason, we defer the construction of the porosity 
change formula beyond the hit-and-stick regime to the fu- 
ture work. We also neglect the rotation of aggregates dur- 
ing their collision for simplicity. Th e rotation could make 
the collisional outcome l ess compact (iRrause & Blumll2004t 
IPaszun & Dominikll2006l) . 

The main goal of this section is to provide a reliable recipe 
for determining the porosity change of aggregates upon gen- 
eral hit-and-stick collisions. As in the previous sections, the 
porosity change due to a single collision are represented by 
the volume of the resultant aggregate, V\+2- Without loss of 
generality, we may rewrite this in the form 

V 1+2 = V 1 + (1 + X)V2, (27) 

where V\ and V 2 (< Vi) are the volumes of the collided ag- 
gregates, and x is a dimensionless factor depending on the 
properties of the aggregates. Since \ vanishes for the ideal 
compact aggregation, we can think of %V 2 as the volume of 
the "voids" newly created in the resultant aggregate. For this 
reason, we may regard \ as the void factor. As we see later, 
the void factor x is constant in the BCCA and BPCA limits, 
suggesting that \ rather than Vi +2 is a "good" variable that 
describes the porosity change due to a general hit-and-stick 
collision. Our task is thus to present an empirical formula for 
the void factor \ m the hit-and-stick regime. 

It is essential to prepare reliable empirical data on which 
the porosity change formula is to be built. For this purpose, 
we have performed numerical experiments on various types of 
collisional aggregation. Although classical BCCA and BPCA 
are useful aggregation models, these are just the two limiting 
cases of general hit-and-stick collisions. In order to bridge the 
gap between the two opposite limits, we introduce a new ag- 
gregation model: quasi-BCCA (QBCCA). We define QBCCA 
as a sequence of ballistic collisions between two clusters with 
a. fixed mass ratio e = M 2 /Mi(< 1). Figure [3] schematically 
illustrates the QBCCA model as well as classical BCCA and 
BPCA. In the following subsections, we describe these aggre- 
gation models in more detail. 

Here we outline how we construct the collision model in 
this section. In Section 4. 1, we introduce how the "volume" of 
a porous aggregate is defined in our model. In Section 4.2, we 
show the results of our Af -body simulations to see how the vol- 
ume of an aggregate evolves in BCCA, BPCA, and QBCCA. 
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(a) BCCA 
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N= 1 



(b) BPCA 
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A/= 1 2 3 
(c) QBCCA (6 = 0.6) 



A/= 1 



Fig. 3. — Schematic description of various types of dust aggregation 
((a): BCCA, (b): BPCA, (c): quasi-BCCA (QBCCA) for e = N 2 /Ni = 
0.6). The filled circles represent particles newly added to the aggre- 
gate at each step. 

In Section 4.3, we synthesize these A^-body results to obtain 
a single empirical formula that approximately determines the 
value of x f° r a general hit-and-stick collision. Finally, in 
Section 4.4, we describe how the aerodynamical and colli- 
sional cross sections of a porous aggregate are determined in 
our model. 

4. 1 . The Characteristic Radius and Volume 

Since a porous aggregate is generally irregular, the concept 
of its "volume," or "radius," is somewhat ambiguous. A useful 
definition for the radius is the gyration radius 



\ 



k=l 



(28) 



where Hk(k = l,2,---,N) are the coordinates of constituent 
monomer and X = N~ l J2k x k * s me coordinate of the cen- 
ter of mass. Anot her definition is the characteristic radius 
(Muk ai etal]ll992h 



a c : 



\ 



3N 



$>*-X)2 



(29) 



Using a c , we can define the characteristic volume 

A<rr 

V c =—a 3 c . (30) 

The characteristic radius has a property that i t reduces to 
the p hysical radius for a homogeneous sphere ( Muk ai et al.l 
1992). In our model, we regard a c and V c as the "radius" 
and "volume" of a porous aggregate 4 . Note that the collision 
model of OST07 uses the area-equivalent volume Va (Equa- 
tion d2TT i) and its associated radius (Equation d26li). 

4 For monomers (N = 1), we simply set a c = uq and V c = Vq. 



Figure [4] shows three samples of A^-body clusters created 
in our numerical experiments together with their "radii" mea- 
sured in different ways 5 . Each of the samples is composed 
of w 10 3 monomers, and hence has the compact (mass- 
equivalent) radius of a* = N l ^ 3 ao « lOarj. The large circles 
in this figure represent the characteristic volume V c of these 
samples. We see that a c well approximates the maximum dis- 
tance from the center of mass to constituent monomers, but 
cia does not. This fact motivates us to define the collisional 
cross section using a c rather than a^ (see Section 4.4.2). 

4.2. Porosity Evolution in Various Types of Aggregation 

Here we summarize our A^-body experiments on three dif- 
ferent types of aggregation (BCCA, BPCA, and QBCCA) and 
derive the volume factor \ for each of the aggregation models. 

4.2.1. BCCA 

BCCA is defined as a sequence of successive collisions 
between identical clusters (see Figure Oa)). In the A^-body 
experiments, we have simulated 10 5 growth sequences of 
BCCA. At each collisional step, we have randomly deter- 
mined the impact parameter as well as the relative orientation 
of colliding aggregates, and followed the ballistic trajectory 
until one of the constituent monomers hit to another. Figure 
shows ten examples of the evolution of a c as a function of 
N. In the figure, we also show the average of a c over 10 5 runs. 
For N > 30, the averaged a c is found to obey a single power 
law 

a c ozN 1/DBCCA , (31) 

where £>bcca is the fractal dimension of BCCA. Our data fit- 
ting shows Dbcca ~ 1 .90, which is consiste nt with the finding 
of previous studies (e.g., Mu kai et alJl!992l) . 

With the power-law relation OTb . it is easy to find the void 
factor x f° r BCCA. Using equations d30l and d3~Tl i. V c ,\+2 is 
written as 



Vr 



N l+2 \ 
N 2 J 



3/.Dbcca 



y e , 2 w 2.99V, 2, 



(32) 



where A^i+2 = N\ +N 2 = 2N 2 . On the other hand, Equation d27l i 
implies V Ct \+ 2 = (2 + xWc, 2 for BCCA. Comparing the two ex- 
pressions, we find 

X « 0.99 = xbcca (33) 
in the power-law limit. Note that Xbcca is independent of V c , i . 

4.2.2. BPCA 

BPCA is a sequence of successive collisions between a 
cluster and a monomer (see Figure Ob)). In the A^-body ex- 
periments, we have simulated 10 3 growth sequences of BPCA 
and obtained the averaged relation between a c and N, which 
is plotted in Figure|3b). For Af > 30, the a c -N relation obeys 
a power law with a fractal dimension w 3.0, and the poros- 
ity P = 1 -V*/V c approaches to P B pca ~ 0.87 4 This result 
is con sistent with the finding of previous work dKozasa et alj 
[HH. 

Using this result, the relation between V c and Af is written 

Vo „ 



V r 



-N. 



1 _ ^BPCA 1 _ ^BPCA 

5 In this paper, the projected area A of a cluster is calculated as the average 
over 15 randomly chosen orientations. 
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QBCCA (e= 




BPCA 




10.1a 
a A =22.9ao 
h a c =32.7ao 



i a 3 =10.1a 

1 a/\=18.5a 

— i a c =19.2ao 

Fig. 4. — Projection of three-dimensional BCCA (left), QBCCA (e = N 2 /N[ = 0.048; center), and BPCA (right) clusters obtained from our numerical 
experiments. Here a,, a A , and a c denote the compact (mass-equivalent), area-equivalent, and characteristic radii, respectively. The gray circles 
represent sphere of radius a c centered on the center of mass. 
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Fig. 5. — (a) Evolution of the characteristic radius a c of BCCA clusters as functions of the number of constituent monomers N. The gray curves 
indicate 10 examples of N-body calculation, and the average over 10 s runs is shown by the black curve, (b) The solid curves show the averages of 
a c (N) for BCCA (top) and BPCA (bottom) limits. The averages for QBCCA models of different mass ratios e =N 2 /N l are indicated by the dashed 
curves. 



Since = Ni + 1, this equation implies 

Vo 



V, 



c.l+2 : 



(35) 

On the other hand, equation ( f27T > implies V c .\+2 = V c ,i +(1 + 
xWo f° r BPCA. Therefore we obtain 



1 



X = 



1 _ ^BPCA 



1 w 6.94 = xbpca, 



(36) 



for the power-law limit. Xbpca is independent of V c ,\ and V c j. 
Equation (l36T l is valid as long as N > 30, or V c \ > (30/1 — 
Pbpca)Vo « 240V . 

4.2.3. QBCCA 

QBCCA is defined as a sequence of successive collisions 
between two clusters with a fixed mass ratio e(< 1). It is clear 
that BCCA corresponds to QBCCA of e = 1 . It is also true that 
BPCA asymptotically approaches to QBCCA of e — > in the 



limit of N\ — > oo. Therefore, this type of aggregation can be 
considered as a model between the BCCA and BPCA limits. 

We here describe the general procedure of QBCCA. Let us 
refer to the larger cluster as the "target," and to the smaller 
one as the "projectile". A target is always chosen to be the 
outcome of the latest collision. A projectile is chosen from the 
outcomes of previous collisions so that the mass ratio between 
the target and the projectile becomes the closest to e. Note that 
this procedure is identical to that of BPCA when A^i < 1 .5/e, 
since the projectile is then always a monomer. In order to 
obtain a "truly" QBCCA cluster, one needs to repeat the above 
procedure until the cluster grows beyond jV w 1/e. 

Figure [3] illustrates an example of QBCCA in the case of 
e = 0.6. The first step is the collision between two monomers, 
as is for any e. The second collision is between the resultant 
dimer (A^i = 2) and a monomer (N2 =1), since N2 = 1 is nearer 
to N\e = 1.2 than N2 = 2. The third collision is between the 
resultant trimer (A^i = 3) and a dimer (A^2 = 2), since N2 = 2 is 
the nearest to A^e = 1.8 among N2 =1,2, and 3. 
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Fig. 6. — Fractal dimension Dq BC ca of QBCCA clusters for different 
mass ratios e = N2/N1 . The value of e = 1 corresponds to the fractal 
dimension Dbcca for BCCA. 

In the A^-body experiments, we have chosen eight values 
of e from the range 0.005 < e < 1.0. For each value of e, 
we have simulated 100-2000 growth sequences and obtained 
an averaged a c -N relation. Figure |5jb) shows two examples 
of the averaged relations (e = 0.05 and 0.005). The averaged 
relations roughly obey a power low 



i c (xN l/DQBCCA(e \ 



(37) 



for N > l/e, where Z)qbcca(c) is the fractal dimension of 
the QBCCA clusters and depends on the mass ratio e. For 
each value of e we measured Dqbcca by fitting the power 
law d37j to the corresponding curve in Figure [5jb). Detailed 
inspection shows that the averaged curves for QBCCA log- 
arithmically oscillate around the power-law fits with a cycle 
AlogAf ~ log(l / e). This is an imprint of the history that clus- 
ters of N < l/e are not constructed by the "truly" QBCCA 
process. We neglect this oscillation in the following modeling 
since the amplitude is small. 

The obtained values of Dqbcca as well as Z)qbcca(I) = 
^bcca are plotted in Figure [6] Remarkably, Z)qbcca(e) dif- 
fers from Dbcca by only 10% even for e ~ 10~ 2 , or V C fl./V Ct i ~ 
(10~ 2 ) 3 ' Dqbcca ~ 10~ 3 . This suggests that an aggregate grow- 
ing by hit-and-stick collisions tends to have a fractal dimen- 
sion close to 2 unless its collision partners are much smaller 
than itself. A simple linear extrapolation of the curve in Fig- 
ure [6] suggests that DqbccaW will approach to 3 if e (or 
V c ^/y c ,\) is set to be as small as 10~ 7 . 

Equation (|37| | implies a relation similar to equation I 



V*. 



1+2 



3/£>qbcca(<0 



= (1 + e) 3 / D QBccA(e) 



(38) 



where we have used N\+2 = (l + e)Ni. We can express this 
equation as a function of V Ci2 /Vc,i instead of e. First, we note 
that equation d37l i also implies 



Yzi. 

V c ,l 



Dqbcca(<0/3 



Since Dqbcca is a function of e, equation 



(39) 

determines e, 



and consequently Dqbcca, as a function of V Cl 2/V Ct i. Using 
equation d39l , equation (1381 is written as 



V, 



c.1+2 : 



1 + 



Ye* 

V C A 



Oqbcca/3 - ] 3/Z)qbcca 



V C A 



(40) 



Thus, the void factor \ f° r QBCCA is given by 



1 + 



'Y* 

:Xqbcca(V c , 2 /V c ,i). 



a/3' 



3/Oq 




(41) 



Since Dqbcca is a function ofV c ^/V c ,i, Xqbcca only depends 
onV C|2 /Ki. 

Figure |7] shows the void factor xqbcca as a function of 
2V c i/(V C A + V c<2 ) = 2/(V cA /V c , 2 + 1). We see that xqbcca 
scales linearly with 2V C 2/(Ve,i + Vc,2)- By data fitting, we ob- 
tain a simple empirical formula 



XQBCCAOWVki) ~ Xbcca- 1.03 In 



V c ,i/V C 2+l 



(42) 



where xbcca is given by equation fl33l l. Note that this for- 
mula includes the BCCA limit since it gives xqbcca = Xbcca 
for V c ,\ = V c ,2- Also note that, although xqbcca increases as 
V c .2 decreases, the void volume xqbcca V 2 decreases with de- 
creasing V c> 2. 

4.3. A Recipe for Porosity Change Due to General 
Hit-and-stick Collisions 

Figure [3b) summarizes xbcca, Xbpca, and xqbcca as a 
function of 2V 2 /(Vi +V 2 ) = 2/(Vi/V 2 + 1). Note that xbpca is 
only plotted in the range V C fl/V c i < 1/240, since the a c -N 
relation for BPCA obeys a power law only when V c j,/V c i < 
1 /240 (see Section 4.2.2). 

Now we combine the above results to construct a recipe of 
the void factor \ f° r general hit-and-stick collisions. In the 
previous subsection, we have seen that \ can be simply writ- 
ten as a function of only V c ^./V c ,\ under a fixed aggregation 
process (i.e., BCCA/QBCCA/BPCA). This can be done be- 
cause N2/N1 is related to V c ,2/V c ,i for each process. This sim- 
plicity is appealing for the modeling of the porosity change. 
For this reason, we assume that the growth history of an ag- 
gregate is well approximated either by successive collisions of 
a fixed volume ratio V c fi/Vc,i (BCCA/QBCCA) or by succes- 
sive collisions of a fixed projectile volume V c 2 (BPCA). This 
assumption allows us to choose \ between xqbcca(Vc, 2 /Vc,i) 
and xbpca- We will discuss the validity of this assumption in 
Section 5.3. Next, we determine the choice so that the value 
of x reduces to Xbcca and Xbpca in the limit of V c , 2 /V c a -> 1 
and 0, respectively. Among the simplest ones, we propose the 
following formula: 

x(Vc,2/V c a) = min{xQBCCA(V c , 2 /V e ,i), Xbpca}, (43) 

where xbpca and xqbcca(V c .2/V c ,i) are given by Equa- 
tions ( f36b and d42l . respectively. This final formula is il- 
lustrated in Figure |7Jb) by the thick grey line. It should be 
noted here that the above choice underestimates the porosity 
increase for QBCCA with very small V c z/V c i. However, as 
seen in Section 5.2, this effect is negligibly small. 

It is wo rth tryin g to compare our hit-and-stick model with 
that of the OST07. Unfortunately, we cannot make a rigorous 
comparison between the two models here since the IOST07I 
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Fig. 7. — (a) Void factor x for QBCCA clusters as a function of 2V c ^/(V c ,i +V Ct2 ). Filled circles indicate the numerical results for various values of 
the mass ratio e. The fitting formula (Eguation (J2j) is shown by the solid line, (b) x for BCCA (Equation £33); filled circle), BPCA (Equation (36l ; 
dashed line), and QBCCA (Equation (42); solid line). The BPCA line is plotted for V c 2 /V c [ < 1/240, for which equation (36} is valid. The gray line 
indicates the void factor formula (43) for general collisions adopted our model. Shown by the dotted curve is the area-equivalent void factor xa 
corresponding to the hit-and-stick model of OST07 (Equation (44)) as a function of 2Va,i/(Va,i +V a ,2>- 



model considers the porosity increase using Va- However, it is 
possibl e to derive the "area-equivalent void factor" \a f° r the 
IOST07I model by equating their hit-and-stick formula d23l and 
V A ,\+2 = V A> \ + (1 +Xa)V a ,2- Neglecting V'add in Equation (|23j 
(which is negligible as long as N\,N2^> 10 2 ), we obtain 



XA,OST07(Va,2/Va,i) = 



Va 



V, 



A. 2 



1 + 



V, 



A. 2 



V 



.4.1 



25cca/3 



1 



1, (44) 



Note that xa,osto7 depends on Va,2/Va,i only. In Figure[7Ib), 
we overplot Xa,osto7 as a function of 2V a ^/(Va^ + Va,2)- We 
see that Xaoswj is approximately equal to x of our model 
in the BCCA limit, but decreases for smaller Va,2/Va,\- As 
shown in Section 5, the OST07 model considerably under- 
estimates the porosity increase of aggregates because of this 
behavior. 

Summarizing Sections 4.2 and 4.3, we have constructed the 
porosity change recipe for general hit-and-stick collisions in 
accordance with numerical experiments on various types of 
aggregation. In the numerical experiments, a new aggregation 
model, which we have referred as QBCCA, has been used to 
fill the gap between the conventional BCCA and BPCA mod- 
els. The porosity change recipe has been given in the form of 
a formula for the "void factor" x (Equation j43l), which even- 
tually determines the volume of a collisionally formed aggre- 
gate, Vi+2, via Equation d27l >. This recipe thus allows us to 
evaluate the porosity change of aggregates upon a single hit- 
and-stick collision. One can calculate the growth and porosity 
evolution of an ensemble of dust aggregates consistently by 
implementing this recipe to the Monte Carlo methods or the 
extended Smoluchowski method presented in Section 2. 

4.4. Cross Section Formulae 

The porosity of aggregates affects the aerodynamical (pro- 
jected) cross section A and the collisional cross section ct co u. 
Below, we describe how we calculate these cross sections in 
our collision model. 



4.4. 1 . Projected (Aerodynamical) Cross Section 

The projected cross section A is one of the most important 
properties of an aggregate since it determines how strongly 
the aggregate is coupled to gas environment. 

It is useful to see how the projected cross section behaves 
differently in BC CA and BPCA. In th e BCCA limit, A is weU 
approximated by (Minato et al. 20061) 



Abcca f 12.57V ' 685 exp(-2.53/N olmo ), N < 16 



A 



0.352/V+0.566/V 



0.862 



7V> 16, 



(45) 



where Ao = na^ is the geometrical cross section of a monomer. 
To see how accurately this formula predicts the actual value 
of A, we take the BCCA cluster in Figure |4] as an example. 
This cluster contains TV = 1024 monomers, and the angle aver- 
age of its projected area is A = 1890a 2 , (the area-equivalent 
radius being = 24.5«o)- The empirical formula, Equa- 
tion d45l >. gives Abcca = 1830a 2 , (a^ = 24.1ao), which agrees 
with the above actual value with a relative error of 3.3%. We 
remark here that the area-equivalent radius a A of a BCCA 
cluster is generally much smaller than its characteristic ra- 
dius a c . For example, the BCCA cluster shown in Figure [4] 
has a c = 61.1ao, which is about three times larger than its 
aA- In fact, the "characteristic" cross section of a BCCA 
cluster, 7rfl 2 «/V 2 / Dbcca Ao«/V l05 A , increases faster than the 
total geometric cross section of constituent monomers, NAq. 
This means that if we naively assumed A = ira 2 c in our colli- 
sion model, we would have a BCCA cluster more and more 
strongly coupled to gas as its growth. In contrast, the formula, 
Equation d45l ), guarantees that A does not increase faster than 
NA . 

In the BPCA limit, A is simply related to a c as 



Abpca ~ ™ c 



(46) 



or equivalently, a a ~ a c . For example, the BPCA cluster 
shown in Figure [4] has the angle -averaged projected cross 
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section of A = 1080oq, or the area-equivalent radius of a a = 
18.5flo. The characteristic radius of this cluster is a c = 19.2ao, 
and thus the relative difference between a a and a c is only 4%. 

Using the above facts, we construct a formula for A of a 
general porous aggregate. Let us take the general formula to 
approximately recover Equations (l45l l and (l46l ) in the BCCA 
and BPCA limits. As one of the simplest ones, we consider 



A(N,a c ) = 



1 



1 



1 



A B cca(A0 na 2 c 7ra~ BCCA (jV) 



(47) 



where Abcca(^V) is the fitting formula defined by Equa- 
tion ( f45b and a C) BCCA(^0 is the averaged characteristic radius 
of BCCA clusters of the same monomer number N. This for- 
mula clearly reduces to Equation ( |43T > in the BCCA limit, and 
also recovers Equations ( f46b in the BPCA limit for large N 
where a? <C A B cca and a c <C fl e .BCCA- The upper panel of 
Figure [8] shows the accuracy of this formula for more gen- 
eral types of clusters. In this figure, the solid curves represent 
the averaged relation between the mass N and the mass-to- 
area ratios N/A obtained for various BCCA/QBCCA models 
(e =1, 0.325, 0.1, 0.05, and 0.01). The dashed curves are ob- 
tained from Equation d4Tb with the averaged a c -N relations 
for our BCCA/QBCCA clusters. For N < 10 6 , the relative 
error between the measured values and the prediction from 
Equation (@7]) is less than 20% for e > 0.05 and less than 30% 
even for e = 0.01. Thus, Equation ( f47T > successfully converts 
a c into A. 

Recently, Pasz un & Dom inik (2009) have proposed another 
simple relation between and the outer radius a ut(~ a c ), 
which reads (see Equation (5) in Paszun & Dominik 2009) 



N 



= 1.21 



^out 

a A 



-0.3 



AT 



-0.33 



(48) 



or A oc (/y L33 a[| u 3 t ) 2 ' 3 ' 3 . They obtained this formula from rela- 
tively small (N < 10 3 ) numerical aggregates. The lower panel 
of Figure [8] examines the accuracy of this formula for larger 
N. The solid curves are the same as those in the upper panel 
of Figure|8] but the dashed curves are now the prediction from 
Equation d48l with q out approximate d by a c . We see that the 



formula of lPaszun & Dom inik ( 2009) underestimates the pro- 
jected area of large (N > 10 3 ) BCCA/QBCCA clusters, espe- 
cially for e > 0.05. 

4.4.2. Collisional Cross Section 

We explained in Section 4. 1 that the characteristic radius a c 
of a porous aggregate well approximates the maximum dis- 
tance from the center of mass to the constituent monomers. 
It will be reasonable to regard two aggregates as collided if 
the impact parameter is smaller than the sum of their charac- 
teristic radii. For this reason we model the collisional cross 
section <7 C0 n for two porous aggregates as 



Ccoii = K(a c \+a c z) 



(49) 



where a c j and a C 2 are the characteristic radii of th e aggre- 
gates. Note that this choice differs from that of OST07, 

Ccoll = 7t(<2a,i +a A ,2) 2 - 

We remark here that Equation (|49l ignores the possibility 
that a collision can be missed even if the impact parameter 
is less than a c i + a e .2- Such collision misses frequently oc- 
cur when colliding aggregates are so fluffy that their fractal 
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Fig. 8.— Mass-to-area ratio N/A of BCCA and QBCCA clusters as a 
function of monomer number N. Upper panel: the solid curves show 
the directly measured values averaged over 10 2 — 10 s samples for var- 
ious values of e (=1, 0.325, 0.1, 0.05, and 0.01 from bottom to top). 
The oscillation seen in the curves for small e is an imprint of the BPCA- 
like growth history for N < 1/e (see Section 4.2.3). The dashed curves 
show the values calculated from our simple formula A(N,a c ) (Equa- 
tion (47}) together with the averaged a c -N relations. The good agree- 
ment between the solid and dashed curves implies that the projected 
cross section is well modeled by Equation (47) as a function of a c 
and N. Lower panel: the solid curves are the same as those in the 
u pper panel, but the dashed curves are obtained from Equation (5) 
of Paszun & Dominik 1 2009, Equation (jjp in the present paper) with 
flout =o c - The formula bv lHaszun & Dominik 12009) underestimates the 
real projected area for very large N. 

dimensions are much lower than two. For ex ample, a micro- 
gravity ex perimen t of lKrause & BlurrJ ((2004) implies that the 
choice of OST07 agrees the actual collisional cross section 
better than ours when the fractal dimension D of aggregates is 
as low as 1 .4. However, the neglect of collision misses does 
not cause a serious overestimation of a co \\ in our model since 
we only deal with aggregates of D > 1.9. 

5. COMPARISON OF THE HIT-AND-STICK COLLISION 
MODELS 

In Section 5, we have constructed a new collision model for 
porous dust aggregates. To compare this model with the clas- 
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TABLE 2 

Summary of hit- and-stick models 



Model 


V 


v l+1 


A 




Compact 


V c 


V c ,l+V c ,2 


A N 2 / 3 


7r(a c .i +a c<2 f 


OST07 


V A 


Equation )23t 


A {V A /V ) 2 / 3 


n(flA,i +a A ,2) 2 


Our model 


Vc 


Equations (27) 


i (43) Equation (47) 


7r(a c ,i+a Cj 2) 2 



sical compact model and the model by OST07, we apply the 
three models to solve simple problems on porous dust coagu- 
lation. 

We consider two types of coagulation problems. In the first 
type, the coagulation is purely driven by Brownian motion. In 
the second type, the coagulation is driven by Brownian mo- 
tion plus relative velocity is of the form Am = gjr/,1 -t/,2|, 
where g is a constant and r/,i and r/,a are the stopping times 
of two colliding aggregates. For example, sedimentation of 
aggregates under a gravitational field g leads to this type 
of relative velocity. Another important example is the rela- 
tive motion driven by turbulence in the stron g coupling limit 
(IWeidensch illing 1984t lOrmel & Cuzzil2 007). where g means 
the typical acceleration of the smallest turbulent eddies. In 
the following, we refer to relative motion leading to the above 
form of Am as differential drift. 

We solve the above coagulation problems using the ex- 
tended Smoluchowski equations ( TTBl l and ( TT6b . Table ^sum- 
marizes the three hit-and-stick models to be compared in this 
section. Here V refers to the volume to be evolved by Equa- 
tion ( TToT i. Note that the characteristic volume V c is identical to 
the area- equivalen t volume Va in the classical compact model. 
Since the OST07 model gives no relation between V c and Va, 
we use the projected area A when we compare the porosity 

evolution of the three models. We set a e ,BCCA(A0 = aoN 1,/DBCCA 
when we use Equation d47b in our model. 

5.1 . Coagulation Driven by Brownian Motion 

The mean relative velocity driven by Brownian motion is 
written as 



Am b = u B0 



1 1 



(50) 



where m B o = ^/clk^T /nntQ is the thermal velocity of a 
monomer. Amb is independent of the projected area A of 
colliding aggregates, so th e area formul a in ou r model does 
not affect the coagulation. Kemp£etaT| (11999b performed a 
full TV-body simulation of Brownian-motion-driven coagula- 
tion and found that the resultant aggregates have fractal di- 
mensions D w 1.8 ± 0.2. We use this fact to examine the va- 
lidity of our collision model. 

In the numerical simulations, we started with a monodis- 
perse ensemble of monomers and followed the evolution of 
n(M) and V(M) using Equations © and (J7J. The control pa- 
rameters are set to Mbd = 80 and S = 0.05. 

Figure [9] shows the normalized mass distribution functions 
N 2 n{N)/no at time t = 10 4 fB0 obtained from the three colli- 
sion models. Here no and fg = no-Ka^uno are the number 
density and growth rate of monomers at the initial moment 
among the models. Note that the mass conservation ensures 
J(N 2 n(N)/n )dlnN =1. We find that the shape of the dis- 
tribution function is insensitive to the choice of the models. 
However, the growth rates of the aggregates are quite different 
among the models. At this time, the mass-weighted average 
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Fig. 9. — Snapshot of the mass distribution function at r = W 4 t B0 for 
Brownian-motion-driven coagulation obtained from different hit-and- 
stick models. The open circles indicate the mass-weighted average 
masses <M),„. 
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Fig. 10. — Snapshot of the p c -M relations at t = 10 4 f s0 for Brownian- 
motion-driven coagulation obtained from the compact model and our 
model. The open circles indicate the average masses M = {M),„. The 
dashed line shows the power-law fit to the density curve of our model, 
and D is the best-fit fractal dimension. Shown by the light gray curve 
is the pc-M relation for our model at earlier time, t = W 2 t BQ . 



mass (M ) m is 10 5 7 mo, 10 6 9 mo, and 10 9 5 mo for the compact, 
OSTOl and our models, respectively (see Section 3.3 for the 
definition of (M) m ). This is essentially due to the difference in 
the porosity evolution among the three models (see Figure fTTI 
below); namely, the higher porosity results in the larger colli- 
sional cross section, and thus the faster growth rate. 

Figure [10] shows the snapshots of the mean internal den- 
sity p c (M) = M/V C (M) at t = 10 4 fB0 for the compact model 
and our model. We find that the density curve for our model 
is well represented by a single power law p c oc M I-3 / D with 
D = 1.96. This means that the aggregates resulting from our 
model have a fracta l dimension of 1.96 , which is consistent 
with the finding of Ke mpf et al.l (1 1 9991) . Thus, our model 
successfully reproduces the porosity evolution of aggregates 
growing in Brownian motion. We also plot the density curve 
for our model at earlier moment, t = 10 2 ?bo- The curve again 
lies on the same fractal line, meaning that the fractal dimen- 
sion of the aggregates does no t change with time. 

The difference between the OST07 model and ours is best 
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Fig. 1 1 . — Snapshot of the normalized mass-to-area ratios N/(A/A Q ) 
at t = 10 4 f B0 for the Brownian-motion-driven coagulation with different 
collision models. The open circles indicate the mass-weighted aver- 
age masses (M),„. The dashed curve shows the mass-to-area ratio in 
the BCCA limit obtained from Equation (45). Shown by the light gray 
curve is the M/A-M relation for our model at earlier time, t = l0 2 t B o- 

illustrated by the evolution of projected area A. Figure QTJ 
shows the normalized mass-to-area ratio N/(A/Aq) for the 
three co llision m odels at t = 10 4 ?bo as a function of N. We see 
that the IOST07I model gives a considerably small projected 
area for growing aggregates co mpared with our model. At this 
moment, the aggregates of the OST07 model with mass (M),„ 
have a projected area approximately an order of magnitude 
smaller than our agg regates o f the same mass. The underes- 
timation of A in the OST07 model is clearly caused by the 
decreasing behavior of Xa,osto7 for Va,i/Va,\ < 1. Thus, the 
realistic modeling of the porosity change between the BCCA 
and BPCA limits is critical. 

5.2. Coagulation Driven by Differential Drift 

We consider aggregates smaller than the mean free path of 
the ambient gas, and adopt the Epstein's law t/ ocM/A. Then, 
the relative velocity induced by differential drift is written as 



Amd = mdo 



N 2 



Ai/Aq A 2 /A 



(51) 



where ujxi = ff r ff) 1& me drift velocity of monomers. 
IWurm & Bluml (0*9981) conducted a laboratory experiments on 
dust coagulation in a turbulent gas environment and showed 
that the resultant aggregates have a fractal dimension of D w 
1.91. 

As in Section 5.1, we start the simulations with a monodis- 
perse ensemble of monomers. However, pure differential 
drift cannot drive the first step of the coagulation because 
Amd vanishes for all pairs of monomers, To avoid this, we 
add Amb as a small perturbation and take the relative veloc- 
ity to be Am = y/ Amq + Am| with u B q = O.Iwdo- The inclu- 
sion of Amb does not affect the dust evolution at a later time 
since Amb decreases as the aggregates grow. The simula- 
tions are performed with the control parameters Afbd = 80 and 
S = 0.01. Time evolution is followed until t = 20?do, where 
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Fig. 12. — Snapshot of the mass spectra at t = 20fo for differential- 
drift-driven coagulation obtained with different collision models. The 
open circles indicate the positions of the mass-weighted average 
masses (M),„. 



fDO = (nonaluDoT 1 . 

Figure [12] compares the mass distributions N 2 n(N)/no for 
the three collision models at t = 20fno- The average mass 
(M) m at this moment is 10 13 8 m , 10 101 m , and 10 8 6 m for 
the compact, IOST071 and our models, respectively. In con- 
trast to the Brownian motion case, our model results in the 
slowest increase in (M) m among the three models. This is 
because our aggregates tends to conserve the mass-to-area ra- 
tio N/A, and thus keep Auo oc A(N/A) small, throughout the 
growth. 

Figure [13] shows the snapshot of the p c -M relation at t = 
20?do for the compact model and our model. We see that the 
density curve for our model behaves differently between the 
low-mass (M < (M) m ) and high-mass (M > (M) m ) sides. At 
the low-mass side, the density curve again obeys a power-law 
relation. The fractal dimension for this s ide is found to be D = 
2.09, which agrees with the finding of IWurm & Bluml dl998l) 
within an error of 10%. The fractal dimension obtained here 
is slightly higher than that in the Brownian motion case. This 
is because the differential drift inhibits equal-sized collisions 
(see Equation ( BP )) and the dominant collision mode shifts to 
lower e. Comparing the obtained fractal dimension D = 2.09 
with Z3qbcca(c) in Figure |6] the mass ratio of the dominant 
collisions is estimated as e w 0.05. We discuss, in more detail, 
the dominant collision mode in the differential drift as well as 
Brownian motion in Section 5.3. 

At the high-mass side, in contrast, the density curve devi- 
ates from the above power law relation and tends to flatten. 
We observe this behavior independently of t, as is found from 
the density curve at an earlier time (f = 10?do) plotted in Fig- 
ure [13] This implies that the growth history of massive ag- 
gregates is qualitatively different from that of smaller ones 
when the coagulation is driven by the differential motion. In 
fact, this high-mass behavior can be explained if we suppose 
that (1) typical-mass aggregates grow by the collision with 
similar-sized ones, while (2) the massive aggregates grow by 
the collision with much smaller, typical-sized (M ~ (M) m ) 
ones. For simplicity, let us divide an ensemble of aggregates 
into two classes, one being a small number of massive ag- 
gregates with mass Mi ^> (M) m , and the other being a large 
population of typical aggregates with mass M 2 ~ (M) m . We 
ignore the collision between massive aggregates and assume 
that the both classes of aggregates grow only by colliding with 
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Fig. 13. — Snapshot of the p-M relations at t = 20t D0 for differential- 
drift-driven coagulation obtained from the compact model and our 
model. The open circles indicate the average masses M= {M) m . The 
dashed line shows the power-law fit to the density curve of our model, 
and D is the best-fit fractal dimension. Shown by the light gray curve 
is the p c -M relation for our model at earlier time, t = lO fan. T he thick 
dotted curves show the analytical prediction by Equation (56}. 



typical ones. In addition, we assume that the both classes of 
aggregates grow at the same rate, i.e., 



1 dM 



1 dM 2 



Mi dt M 2 dt 



1 



'grow 



(52) 



where f * w is the growth rate independent of M 6 . Under these 
assumptions, the increase in the volume V 2 of a typical-mass 
aggregate is given by the rate of collisions with other typical- 
mass aggregates, ^,22 = (dM 2 /dt)/M 2 = t^ ov/ , times the vol- 
ume increase per collision, [1 + xOV^z)]^ i- e -, 



dV 2 _ 1+ Xbcca 
dt 



V 2 , 



(53) 



'grow 



where we have used that x(^2/V 2 ) = Xbcca- Meanwhile, the 
increase in the volume V\ of a high-mass aggregate is given by 
the rate of collisions with typical-mass aggregates, ^ C oii,i2 = 
(dM\/dt)/M 2 = {M\/M 2 )t~} ow , times the volume increase per 
collision, [l + x(V2/Vi)]V2, i-e., 



dVi 
dt 



Mi l+x(V 2 /Vi) T7 Pi l+x(V 2 /Vi) 



M 2 



-V 2 



P2 



Vi, (54) 



where we have used that p = M/V. Using Equations ( T52l - 
(|54t . the time variation of p\ / p 2 is written as 



dln(pi/p 2 )_ 1 dM { 1 dVi ( 1 dM 2 
Vy~dT~ \W 2 ~dT 



dt 



Mi dt 
1 dV x 



1 dV 2 \ 

V 2 ~d7J 



_l_dv2 

V\ dt V 2 dt 



= J-(i+Xbcca--[1 + x(V2/Vi)]}. (55) 

'grow I, Pi ) 

The solution to this equation asymptotically satisfies 

6 This assumption works well when the coagulation is driven by the dif- 
ferential drift. In this case, dM/dt = Pjcf^w Ahd scales with (cr co u/A)M, so 
f™ w scales with <r co u/A. However, ct co u is roughly proportional to A, and 
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Fig. 14. — Snapshot of the normalized mass-to-area ratios N/(A/A ) 
at t = 20t m for the differential-drift-driven coagulation with different col- 
lision models. The open circles indicate the mass-weighted average 
masses (M),„. The dashed curve shows the mass-to-area ratio in the 
BCCA limit obtained from Equation (45). The dot-dashed curve is the 
result obtained from our model with setting \ = Xqbcca for all V 2 /Vi 
instead of Equation (43). Shown by the light gray curve is the M/A-M 
relation for our model at earlier time, t = U)t DQ . 



d(p\/p 2 )/dt = 0, or 

1+ Xbcca 



Pi = 



l+xOVVi) 



P2 = 



1 + Xbcca 



1 + X 



M 2 /Mj 
P2/P1 



pi- 



rn 



therefore t ' w is approximately constant. 



For fixed M 2 and p 2 = p(M 2 ), this final equation implicitly 
determines pi as a function of M\. In Figure [13] we over- 
plot the solutions p\ to Equation d56l > at t = IQtm and 20foo 
as a function of M\ with M 2 = (M) m and p 2 = p((M) m ). We 
find that these solutions reproduce the flattening of the density 
profiles. The flattening of the density curve is intuitively rea- 
sonable, since the porosity of an aggregate is generally kept 
small when it only accretes much smaller aggregates. 

Figure [14] compares the normalized mass-to-area ratios 
N/(A/Aq) for the two collision models at t = 20foo as a func- 
tion of M. The mass-to-area ratio is a key property of ag- 
gregates growing by differential drift since the terminal ve- 
locity is proportional to it. Again, the M/A-M relation for 
our model exhibits different behavior between the low-mass 
(M < (M) m ) and high-mass (M > (M) m ) sides: nearly flat in 
the low-mass side and increasing with mass in the high-mass 
side. This difference clearly reflects the p-M relation seen in 
FigureQj] Moreover, detailed inspection shows that M/A for 
the lOST07l model also steepens in the high-mass end. Thus, 
the steepening of the M/A curve is not peculiar to our model. 

One might suspect that the increase in the mass-to-area 
ratio is caused by the cutoff in the void factor x (Equa- 
tion d43b ) at small V 2 /V\. Actually, the complete flatten- 
ing as seen in Figure [13] is achieved because we have cho- 
sen x = Xbpca in the small V 2 /V\ limit (see Section 4.3). 
This can be confirmed from Equation d56T ). which implies 
Pi ->■ P2(1+Xbcca)/(1 + Xbpca) ~ 0.25p 2 independently of 
Vi. However, flattening still occurs even if we choose \ to 
be Xqbpca(V 2 /Vi) for all V 2 /V\. In fact, Equation (l56b then 
implies p\ oc l/log(Vi/V2) in the small V 2 /V\ limit, which 
decreases only logarithmically with increasing V\ . For con- 
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firmation, we have solved the same coagulation problem by 
removing the cutoff in \ and assuming \ - Xqbcca(V2/Vi) 
for all V2/V1. The result is overplotted in Figure [14] We see 
that M/A at the high-mass end is mostly unchanged even if 
the cutoff in x is removed. Thus, the increase of the mass- 
to-area ratio seems to be a robust feature for the coagulation 
driven by the differential drift. 

The above finding has a potentially important implication 
for the growth of dust aggregates in protoplanetary disks. One 
of the authors has recently examined the effect of d ust charg- 
ing on its collisional growth in protoplanetary disks (Okuzumi 
2009). Aggregates in a disk usually charge up negatively on 
average, and their collisions require a sufficiently high col- 
lision velocity to overcome the electrostatic barrier. Strik- 
ingly, it turned out that the electrostatic repulsion strongly 
suppresses the growth of fractal (D w 2) aggregates. This 
is essentially because a fractal aggregate of D w 2 keeps a 
small mass-to-area ratio and thus a low sedimentation veloc- 
ity. However, our calculation shows that high-mass aggre- 
gates in an ensemble have a large mass-to-area ratio and a 
high sedimentation velocity compared with low-mass fractal 
aggregates. This implies that the high-mass aggregates may 
be allowed to continue growing even if the others are not al- 
lowed. In a forthcoming paper, we will investigate this issue 
using the numerical tools presented in this study. 

5.3. Validity and Possible Limitation of Our Porosity 
Evolution Formula 

In the modeling of the porosity evolution, we have as- 
sumed that the growth of an aggregates is well approximated 
as BCCA/QBCCA or BPCA. Here, we examine whether this 
assumption holds in the above coagulation problems. 

Let us consider a quantity 

C Ni (e) = — x , (57) 

/ t'K{e'Ni;Nx)n{e'Nx)de 
Jo 

where K(eNi;N\)n(eN\) is the rate of collisions between ag- 
gregates with mass eN\ and iVi . The growth rate t~* ow = 
(dN\ I dt) /N\ of an aggregate with mass N\ is written as f ~ow = 

Jo' e'K(e'Ni;Ni )n{e'N] )d(e'Ni). Thus, C Nl (e)de measures how 
collisions of mass ratios between e and e+de contribute to the 
growth of an aggregate of mass 2Vi(> eN\). Figure [TBI plots 
C(if\ (e) for Brownian motion and differential drift at differ- 
ent times t. For both the cases, the functional form of C(/v) m (e) 
is nearly independent of t and has a peak at e = 1 (Brownian 
motion case) and e ^ 0.1 (differential drift case). The dis- 
tribution of C/jy\ (e) has an order-of-magnitude spread in e- 
space. However, this spread can be regarded as narrow since 
an order-of-magnitude fluctuation in e only causes ~ 10% 
fluctuation in Dqbcca (see Figure [6J. Therefore, the growth 
of aggregates is indeed well approximated as BCCA/QBCCA 
when the coagulation is driven by Brownian motion and dif- 
ferential drift. 

On the other hand, the assumption may become less good 
if the collision partners of an aggregate have widely spread 
mass ratios and if they contribute equally to its growth. Such 
a situation may be realized when, for example, the collision 
partners obey a very broad and flat size distribution. 

6. SUMMARY 




mass ratio e = N/(N) m 

Fig. 15. — Normalized, mass-weighted collision rates C^ m (N) be- 
tween aggregates of masses (N) m and N = e{N),„ as a function of 
e =N/(N) m (thick solid: Brownian motion, t = l0 4 t B0 ; thin solid: Brow- 
nian motion, t = 10 2 f B0 ; thick dashed: differential drift, t = 20too\ thin 
dashed: differential drift, t = Wt D0 ). Note that the two solid curves are 
indistinguishable. 

In this study, we have presented new numerical tools for 
studying the coagulation and porosity evolution of dust ag- 
gregates. Our findings are summarized as follows. 

1. We have presented a new numerical method for sim- 
ulating the coagulation and porosity evolution of dust 
aggregates as an extension of the conventional Smolu- 
chowski method (Section 2). The new method treats 
the averaged volume of aggregates of the same mass 
as dynamical, and follows its evolution as well as the 
evolution of the mass distribution function consistently 
(Equations ( TT31 > and ([T6ll). This method enables to treat 
the 2D (i.e., mass and porosity being dynamical) co- 
agulation problems in a very efficient way. We have 
confirmed that our method well reproduces the results 
of previous full 2D Monte Carlo simulations with much 
fewer computational expense (Section 3). 

2. We have presented a new collision model based on our 
numerical experiments on aggregates collisions (Sec- 
tion 4). A collision model refers to a set of definitions 
and formulae for the properties of porous aggregates, 
including a recipe for the porosity change upon a col- 
lision. As the first step, we have focused on "hit-and- 
stick" collisions, i.e., collisions involving no restructur- 
ing or fragmentation. In the numerical experiments, we 
have first considered successive collisions between ag- 
gregates of a fixed mass ratio, which we have referred 
as QBCCA. This allows to fill the gap between the clas- 
sical BCCA and BPCA. Using the results of the JV-body 
experiments on BCCA, BPCA, and QBCCA, we have 
constructed a recipe for the porosity change due to a 
general hit-and-stick collision (Equation d43l ) as well 
as formulae for the aerodynamical and collisional cross 
sections. In a forthcoming paper, we will include the ef- 
fects of collisional compression and fragmentation into 
our collision model. 

3. To clarify the validity of our collision model and the 
difference from previous models, we have performed 
a couple of simple simulations on porous dust coagu- 
lation with the extended Smoluchowski method (Sec- 
tion 5). We considered two cases where the coagu- 
lation is driven by Brownian motion and differential 
drift, respectively. For both cases the simulations us- 
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ing our model result in fractal aggregates of D w 2, 
which is consistent with the findin gs of previous full 
Af-body and laborato ry experiments ( Kempf et al.ll999t 
IWurm & Blurrifl998l) . By cont rast, the simulations us- 
ing the hit-and-stick model of OST07 result in much 
less porous aggregates. This is because this model un- 
derestimates the porosity increase upon unequal-sized 
collisions. 

4. We have discovered that, when the coagulation is driven 
by differential drift, aggregates at the high-mass end 
of the mass distribution obey a flat density-mass rela- 
tion, in marked contrast to aggregates of lower masses, 
which obey a fractal (i.e., power-law) density-mass re- 
lation (Section 5.2). This is explained by the difference 
in the dominant growth modes: typical aggregates grow 



by colliding with similar ones, while the most mas- 
sive aggregates grow by colliding with much smaller 
ones. As a consequence of this tendency, the massive 
aggregates can have a drift velocity (oc mass-to-area ra- 
tio) considerably higher than smaller ones. This find- 
ing may be crucially important in protoplanetary disks 
where the growth of slowly moving aggregates is sup- 
pressed due to the negative charging. 



We thank Andras Zsom for providing us with useful infor- 
mation on his Monte Carlo simulations. We also thank the 
referee, Chris Ormel, for fruitful comments leading to signif- 
icant improvements in the manuscript. 
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